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ABSTRACT 

We have identified over 2000 well resolved cluster halos, and also their associated bound subhalos, 
from the output of a 1024 3 particle cosmological N-body simulation (of box size 320/i _1 Mpc and soft- 
ening length 3.2/i -1 kpc). This has allowed us to measure halo quantities in a statistically meaningful 
way, and for the first time analyse their distribution for a large and well resolved sample. We charac- 
terize each halo in terms of its morphology, concentration, spin, circular velocity and the fraction of 
their mass in substructure. We also identify those halos that have not yet reached a state of dynamical 
equilibrium using the virial theorem with an additional correction to account for the surface pressure 
at the boundary. These amount to 3.4% of our initial sample. For the virialized halos, we find a 
median of 5.6% of halo mass is contained within substructure, with the distribution ranging between 
no identified subhalos to 65%. The fraction of mass in substructure increases with halo mass with 
logarithmic slope of 0.44 ± 0.06. Halos tend to have a prolate morphology, becoming more so with 
increasing mass. Subhalos have a greater orbital angular momentum per unit mass than their host 
halo. Furthermore, their orbital angular momentum is typically well aligned with that of their host. 
Overall, we find that dimcnsionless properties of dark matter halos do depend on their mass, thereby 
demonstrating a lack of self-similarity. 

Subject headings: cosmology: dark matter — galaxies: clusters: general — methods: N-body simula- 
tions 



1. INTRODUCTION 

N-body dark matter simulations have proven to be one 
of our most valuable tools for testing the predictions of 
the standard model of hierarchical structure formation. 
Indeed, it is largely through studying the structures that 
form in such simulations that we are able to investigate 
proposed models that describe the relative mass distri- 
bution of halos and the subhalos within them, as well 
as their physical characteristics such as spin, shape, and 
concentration. 

Until the end of the last decade, it was not possi- 
ble reach the required mass resolution to enable the 
identific ation of galaxy mass halos as substructures in 
cluste rs llWhitdll 976Uva,ri Ka,mnenllT995l ISummers et aTl 
119951: iMoore et all Il996|) . Commonly known as the 
overmerging problem, this was mainly due to the limited 
mass and force resolution of the simulations used. The 
major causes of this problem were premature tidal 
disruption due to inadequate force resolution and two- 
particle e vaporation for halo s with a small number of 
particles (Klvoin et al. 1999). However, rapid advances 
in parallel computing, through both the improvement of 
hardware and the development of fast and efficient par- 
allel algorithms, has enabled us to achieve the numerical 
resolution required to overcome the numer i cal problems 
iGhigna et all Il998at iKlvpin et all 119991 IMoore et all 
1999t lOkamoto fc Habel 119991: iGhigna et alJ 120001: 
Sorinerel et alJ 120011 [Bode et al ] 120011 l Kray1;so v et al.l 
2004aI lDeLucia et al.ll2004t ISprineel et al.ll2005fl 

There are two preferred approaches for investigating 
the physical properties of simulated dark matter ha- 
los. In order to achieve high numerical precision — 



thus allowing the accurate identification of low mass 
subhalos and measurement of the properties in the in- 
nermost regions — many previou s studies have a dopte d 
a resimulation techn ique (e.g. INavarro et alJ |l995), 
Navarro et alJ (I1996J), Thereafter referred t o as NFW), 
Gao et all <l2005l) INavarro et all (|2004fK IMoore et all 
l|1999|) . lReed et all <|2005a|) L This involves extracting a 
halo from a cosmological simulation and resimulating it 
at much higher resolution, approximating the tidal field 
of the surrounding structures. By picking a 'nice' (or re- 
laxed) halo in the original simulation, it can furthermore 
be ensured that the final hal o will be in dynamic eq uilib- 
rium. Using this technique, IDiemand et al.l l|2005|) com- 
pleted the first billion particle halo simulation, achieving 
unprecedented resolution and probing the halo density 
profile to within 0.1% of the halo virial radius. 

However, there are two downsides to this approach. 
Firstly, resimulating a small number of halos prevents an 
accurate determination of the distribution their physical 
properties. Previous studies have typically resimulated 
several halos in an attempt to ascertain if geometric and 
dynamical properties differ with mass. There remains 
a significant degree of uncertainty in the results of such 
studies merely because of the small number statistics in- 
volved. The resimulated halos do not present an accurate 
picture of the entire range of physical characteristics that 
may be adopted by dark matter halos when one chooses 
only relaxed halos for higher resolution. 

The second approach that one can take is to uni- 
formly simulate a cosmological volume at the high- 
est resolution possible (given the available resources) 
and pick out a large number of halos. Although the 
ability to achieve higher possible mass resolution is 
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lost, this method results in a statistically significant 
sample of halos. A number of previo u s studies (e.g. 
IDavis et all (fl985h. IWarren et alJ jl99^) . ICole fc Lacevl 
(1996), Bullock ct al. ( 2001,3)) employed this approach 
to analyse the distribution of the integral properties 
of a large sample of halos, including their virial mass, 
spin and overall distribution over the volume simu- 
lated; or they used the dark matter distributions to esti- 
mate the lar gje angle gravitational tens in g of background 
sources (e.g. iWambsganss et al.ll2004h . iKravtsov et al.l 
l)2004aj) performed simulations with high enough reso- 
lution to identify substructures and even subsubstruc- 
tures within the cluster halos extracted from the sim- 
ulation box. They analyze the halo occupation distri- 
bution for the host halos and the two-point correlation 
function of the galaxy size halos and subhalos, finding 
the results to be in good agreement wi th those from 
the Sl oan Digital Sky Survey. Recently, iSpringel et al.l 
(2005) presented initial results from the "Millenium" 
simulation - the highest resolulation cosmological N- 
body simulation to date - containing 2160 3 particles 
in a cube of side length 500/i _1 Mpc. This simulation 
was performed using the new parallel GADGET2 code 
(Springel 20051) which adopts a variant of the Tree- 
Particle-Mesh (hereafter, TPM) technique f or calculat- 
ing the gravitational forces o n each particle ( Bod e et all 
200lt iBode fc Ostrikerll200l . 

This paper takes the second approach, presenting an 
analysis of a 1024 3 particle simulation with box size 
320 h~ x Mpc evolved using the TPM algorithm of 
IBode fc Ostrikerl 1)20031) and containing more than 57,000 
structures at the last timestep. Of these we have anal- 
ysed the most massive 2200 cluster halos containing at 
least 10,000 particles each, thus achieving unprecedented 
resolution in such a large sample of halos. The clusters 
are first identified using a friends-of-friends algorithm; 
we the n employ the geometrically based Denmax rou- 
tine of iBertschinger fc Gelbl l)1991fl to identify the sub- 
halos in each cluster and build family trees. Finally, all 
energetically unbound particles are removed. The algo- 
rit hm we employ up to this stage is described in detail 
in lWeller etaDil200l . 

We first present a procedure to identify the halos that 
are still in the processes of formation and have not yet 
achieved a dynamically relaxed, virialized state. The 
purpose of this is to enable an independent analysis of the 
physical characteristics of those halos that are dynami- 
cally relaxed. A catalogue of the main physical proper- 
ties of the halos is then compiled, including morphology, 
spin, fraction of mass in substructure, the subhalo mass 
function, concentration and maximum circular velocity. 
We analyse the distributions of each quantity, which are 
well defined due to the size of our sample, the depen- 
dence, if any, on halo mass, and how halo substructure 
influences overall cluster properties. The key questions 
we ask therefore, are 

• What is the effect of the dynamical state of the 
halo on the physical properties that we measure? 

• What are range the of values that the physical 
properties of the halos take (in the mass range we 
investigate)? 

• Are cluster halos self-similar, i.e. how do the prop- 



erties of halos vary with halo mass? 

• What is the effect of substructure on the overall 
characteristics of a cluster? 

In Section 2 we describe in detail the sample of ha- 
los that we analyse in this paper and how unvirialized 
clusters can be identified. In Section 3 we discuss the 
physical properties that we measure and how we do so. 
In Sections 4, 5, and 6 we present the main results of 
our analysis — the distribution of halo properties, how 
they depend on halo mass and the impact of substruc- 
ture on halo properties. Finally, in Section 7 we discuss 
the results and state our conclusions. 

2. DEFINITION OF SAMPLE 
2.1. The Identification of Halos and Subhalos 

There are many possible algorithms for defining and 
identifying halos and subhalos, and there is no gener- 
ally agreed upon best method adopted by current inves- 
tigations. Some methods are essentially geometrical and 
others aim at finding dynamically coherent structures; 
all contain both dimensionless and dimensional parame- 
ters. In general, there are three levels of refinement that 
one can adopt. The most basic is to use a geometrical 
routine such as the friends-of-friends llHuchra fe Gellerl 
119821 IDavis et al]IT985l ILa,cev fc Cjoldll 99411 or Denm^ 
(IBertschinger fc Gelbl 119911 iGelb fc Bertschingerl Il99l 
lEisenstein fc Hutlll998() algorithms. These use only in- 
stantaneous particle positions to group together nearby 
particles (defined by the FOF 'linking length' or the Den- 
max 'smoothing length') into localized structures. Nei- 
ther method performs any type of dynamical analysis to 
check whether these structures are bound. 

The next level of refinement is to follow a geometri- 
cal identification of subhalos with a procedure for itera- 
tively removing the unbound particles with the greatest 
total energy until only bound particles remain. Exam- 
ples of routines th at use t his approach i nclud e SKID 
ilStadel et al]|!997l) . BDM llKlvpin et al] Il99l. MHF 
\Gill et al.ll2004aft and SUBFIND ijSnringel et alMM!^ . 
In each of these, the unbinding is performed assuming 
the subhalo is completely isolated, i.e. only the bound 
particles in the halo at each iteration are taken into con- 
sideration when calculating the potential energy of each 
particle. Not taken into account in this energy calcu- 
lation are the (previously identified) unbound particles 
located spatially within the subhalo, nor the disruptive 
effect of the tidal forces from the particles surrounding 
it. In iShaw et al] l)2006j) we investigate the impact of 
allowing for these forces, finding there to be little differ- 
ence in the subhalo populations measured; the increase 
in the binding energy of a subhalo from including all the 
particles located within it is almost entirely balanced by 
the losses due to the external forces. 

In this paper we follow previous authors and imple- 
ment only the first refinement — identifying energeti- 
cally bound structures — to an initial geometrical selec- 
tion of structures. The a lgorithm is described in de- 
tail by iWeller et all l)2005|) and is applied to a ACDM 
N-body simulation of 1024 3 particles with box size 
320fe~ 1 Mpc and a splin e kernel force softening length 
l|Hernauist fc Katzil989l) of e = 3.2ft- 1 kpc. This simula- 
tion was evolved to 2=0.05 using the Tree-Particle-Mesh 
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(TPM) algorithm ijBode fc Ostrikerll2003|) . The cosmo- 
logical parameters used include fi m = 0.3, A = 0.7, 
and <78 =0.95; outputs from this run have previously 
been used to make predictions concerning stro ng lensing 
ijWambseanss et al.l2004UHennawi et al.l2005|) . We wish 
to ensure that all the halos are well resolved and that the 
overmerging problem is not in evidence. Thus we discard 
all halos with less than 10,000 particles. The minimum 
halo mass allowed is therefor e w 3 x 10 13 / i, ~ 1 M|7) . 

The method described in iWeller et alJ l|2005l) starts 
with the hierarchical identification of structures: firstly 
the large clusters are identified through a Friends-of- 
Friends (FOF) routine with a linking length of b = 
0.2n -1 / 3 , where n" 1 / 3, is the me an inter-particle sep- 
aratio n. The Denmax routine of iBertschinger fc Gelbl 
(1991) is then run on each FOF halo with a high res- 
olution smoothing length of 5e in order to identify the 
substructures within each cluster halo. A family tree is 
then constructed by hierarchically associating the small- 
est mass subhalos with their lowest mass 'ancestor', so 
that each subhalo has only a single immediate parent. 
Those that consist of less than 30 particles are dissolved 
into their immediate parent. 

Up until this point the analysis is purely geometrical. 
Next, energetically unbound particles are discarded, first 
by removing the velocity outliers, and then by iteratively 
identifying those particles with a total energy greater 
than zero (in the center of mass frame of the halo) and 
removing the most energetic. If at any stage the mass 
of a subhalo drops below 30 particles, we dissolve it into 
its immediate parent. We then check to see if pairs of 
the immediate daughters of the parent cluster halo are 
bound. If this is the case, they form a hyperstructure: 
the less massive of the two becomes a daughter of the 
more massive structure. Finally, we remove all subhalos 
that are not bound to the biggest structure — the mother 
halo — and any particles that are not bound to the entire 
cluster. As described in IWeller et al.l l)2005|) . this entire 
procedure is both stable and largely independent of ar- 
bitrary parameters choice. Henceforth, we refer to the 
most massive structure in each halo as the 'mother' or 
'host' halo. Substructure particles are those that are also 
bound to any other smaller structure within the halo. We 
use 'cluster' or 'the halo 1 to refer to the entire assembly. 

In this study, we also incorporate an additional step: 
the tagging of 'irregular' halos in our sample. The pur- 
pose of this is to allow a separate analysis of those halos 
that are still undergoing significant changes through their 
involvement in a major-merger. This is described in the 
following Section. 

2.2. Virilization Criterion — Identifying Irregular 

Halos 

One of the main aims of this study is to compare the 
physical properties of the virialized halos in our sample to 
those that are still in the process of forming. To achieve 
this we sort through the halo catalogue to isolate and tag 
two types of halo: 

• those with mass below a certain threshold, to en- 
sure that all the halos are well resolved and that 
the overmerging problem is not in evidence. These 
are eliminated from the sample. 



• halos that we deem to have not yet reached a state 
of dynamical equilibrium. These are tagged so that 
we can analyse their physical characteristics sepa- 
rately. 

To satisfy the first condition, we discard all halos with 
less than 10,000 particles, or w 3 x lO 13 ^^/©. We use 
the virial mass as our measure of halo mass. For a ACDM 
cosmology, it is conventional to define the virial mass 
M v ; r and radius i? v ir as M v ; r = |7ri? 3 ir A c (z)p c (z), where 
p c is the critical density of the universe, and the m ean 
over-density A c = 178fi m (z) 45 i|Lahav et al.lll991fl . In 
order to calculate the virial mass of each halo, we start at 
its density maximum and proceed outward until we reach 
the virial radius, within which the mean over-density is 
A c . Further, we include substructures whose center of 
mass is within the virial radius. Hence we, define a par- 
ticle as being part of the halo if it is bound to the halo 
and lies within the virial radius, or is part of a bound 
subhalo that has its center of mass within the virial ra- 
dius. When quoting halo masses, we do so relative to 
the characteristic mass scale, = 8.0 x 1O 14 /i -1 M0, 
obtained by fitting a Schechter function to the mass dis- 
tribution of all halos with greater than 5000 particles in 
the simulation box (see Fig ure 2 in IWeller et all <|2005)). 
At this stage, having removed the low mass halos, we 
have an overall sample of 2159 halos. 

For the second requirement, we use the virial theorem 
to determine which halos are not dynamically relaxed. 
The virial theorem for a system in equilibrium is usu- 
ally stated as 2To + Wq — 0, where To is the total ki- 
netic energy of the halo, and Wq the total potential en- 
ergy. However, this assumes a system with finite phase 
space volume, and also that each halo is in complete iso- 
lation and that all mass associated with the halo has 
been considered. In selecting a halo, a cut is made at 
some outer boundary, in this case corresponding to the 
isodensity contour at which the spherically a veraged den- 
sity i s A c « 100 times the critical density ijLahav et alJ 
1991). Any particles bound to the halo outside of the 
virial radius are by the standard convention not included 
when we calculate the total gravitational and kinetic en- 
ergy for each halo. However, these particles do make 
a significant overall contribution to the pressure at the 
halo boundary, and so we must incorporate an additional 
term into the virial theorem to account for them. This 
correction manif ests itself as a form o f surface pressure 
at the boundary ijChand rasckhar 196l]) of the halo, with 
the energy content 

E s = J P s (r)r-dS . (1) 

The full and exact version of the virial theorem for a 
self-gravitating system is thus 

l§=^o + Wo-E s , (2) 

where I is the moment of inertia of the halo. We therefore 
define as our measure of 'virialization' the following: 



where j3 — > as J 
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In practice, we compute P s in the following manner: 
First, we rank order all the particles in the halo by radius 
and select the outermost 20%. We label the radius of the 
innermost particle in this shell as Rq.s, the outermost as 
i? v ir and the median as Ro.g- Using the ideal gas law, we 
can approximate the surface pressure term, 



Ps = (1/3) 



E<(' 



V 



(4) 



summing over all particles that lie between R os and R v 
V is the volume occupied by the outermost 20%, 



V 



4nRl 8 



(5) 
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hence Equation n can be approximated by 

E s » 4nR 3 9 P s . (6) 

Figure n compares the virialization distribution of 
all the halos in our sample, omitting and including the 
surface pressure term respectively. At very low redshift, 
it is to be expected that the majority of the dark mat- 
ter clusters in a simulation will be dynamically relaxed. 
By comparing these two distributions it is clear that our 
inclusion of the surface pressure term in the virial the- 
orem is necessary — the distribution in the lower panel 
now peaks very close to zero. The tail of the distribution 
to negative values of (3 is largely produced by systems 
which are infalling and not yet virialized. 

However, the virial equilibrium does not imply that 
d 2 I/dt 2 — instantaneously at all times, but rather that 
< d 2 I/dt 2 >= 0, time-averaged over a period that is long 
compared to the local dynamical timescale. Therefore, 
we expect a roughly symmetric distribution of (3 ± op 
around zero due to those halos that are oscillating about 
the virial equilibrium position. It is difficult to know 
what exact value of op to expect as it depends on the 
state of subdivision of the system. Thus, as there are also 
no halos with (3 > 0.2, we have picked the value of the 
cutoff, op = 0.2, so as to leave the remaining distribution 
of (3 roughly symmetric about zero, as indicated by the 
arrow in the lower panel of Figure^] Although this cutoff 
value is essentially determined arbitrarily, we feel that it 
sufficiently removes the particularly unvirialized halos in 
the tail of the distribution. Applying this criterion to our 
sample results in 74 (3.4%) of our halos being tagged as 
'unvirialised'. In Section 01 we demonstrate that these 
halos have significantly different physical properties to 
those that adhere to our criterium. Reducing the value 
of op dilutes this difference, whereas increasing it does 
not significantly change our results. 

We should also be able to apply Eq. |3] at any radii 
in order to determine whether the mass within is in dy- 
namically relaxed. Figure shows the variation of (3 
with radius for three separate halos, including (circles) 
and omitting (crosses) the surface pressure correction. 
To create this plot, we have calculated the total kinetic 
potential energy for all particles within a certain fraction 
of the virial radius, r/i? v j r , for spheres of successively in- 
creasing radius r out to i? V ir- We calculate the surface 
pressure at each r by applying Equation 0] to the out- 
ermost 20% of particles within that radius. The three 
halos have been selected according to their values of (3 at 
i?vi r : zero (top), 0.1 (middle) and -0.21 (bottom panel), 
using Equation |3 Thus the latter would be identified as 
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Fig. 1. — Histograms of the distribution of /3 for the halos in our 
sample having omitted the pressure term E a from Equation [3] (top 
panel). The distribution peaks at a value just greater than -0.1, 
rather than at as would be expected for virialized halos. Bottom 
panel — same as top but having included the pressure term whilst 
calculating /3 for each halo. The distribution now peaks much closer 
to zero. Also indicated (by the dashed line) is the value at which 
we take our lower limit of /3. Halos with /3 < —0.2 (shaded) are 
removed from our sample for the following analysis. 



being unvirialized, according to our criterium. Density 
plots of the particle distributions for each of these halos 
can be viewed in Figures 

It is clear that the value of 2T/W + 1 (i.e. j3 omitting 
the surface pressure) is negative at all radii for all three 
halos. However, if we had been able to continue the cal- 
culations beyond the virial radius for the halos in the top 
two panels in Figure [3 it is likely that bo th would con- 
verge to zero. We note that this is found bv lCole fc Lacevl 
( 1996) who show in a similar plot that the ratio 2T/| W| 
converges to unity for their halos at large radii (see their 
Figure 8). By neglecting to include the particles out- 
side of the virial radius we are omitting an important 
term in the virial theorem required to demonstrate that 
these halos are in fact in dynamic equilibrium. For these 
two halos, the surface pressure corrections demonstrate 
that the correction achieved through accounting for the 
finite pressure and density at each radii result in values 
of (3 close to zero throughout. These halos appear to 
be uniformly in dynamical equilibrium, at all radii. The 
projected density plots of their particle distributions in 
Figure (labelled a and b) suggest that they are both 
relaxed clusters - a massive central core with smaller or- 
biting substructures. 

The radial distribution of (3 for the third halo (bot- 
tom panel) — also the most massive of three — is con- 
siderably more negative at each radii than the other two 
halos. This halo has not yet reached dynamic equilib- 
rium; even with the surface pressure correction, which is 
small compared to the halo kinetic energy, (3 still remains 
much less than zero at all radii and appears to converge 
to a constant value at r/R v - u = 0.7. Plot c in Figure 
|2| demonstrates that this halo is clearly a complicated 
structure. It appears to have several clumps within its 
core, suggesting that it has recently undergone a major 
merger. 

Figures 0] and show how j3 varies with halo mass and 
with the fraction f s of halo mass in substructure. Both 
of these plots also include lines representing the median 
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Fig. 2. — Virialization parameter j3 (Eq. as a function of radius 
for three example halos, excluding (crosses) or including (circles) 
the surface pressure term E s . 
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values in mass and f s bins. Only halos with j3 > —0.2 
have been used in calculating the trend line. From Fig- 
ure 01 it appears that there is a weak trend for (3 to 
decrease (below zero) as halo mass increases. In accor- 
dance wit h the current model of hierarchical structure 
formation ijLacev fc Coldll993|) . larger and more massive 
objects have formed more recently through the merger of 
smaller structures. Hence, the most massive dark matter 
clusters are still in the process of dynamical relaxation 
and their constituent particles have an excess of kinetic 
energy, which we would expect to dissipate as the cluster 
settles into dynamic equilibrium. Halos with larger val- 
ues of M V i r /M* are, statisically, younger than those with 
smaller values of this quantity and thus are, on average, 
less 'virialised'. The same explanation applies to Figure 
— which shows (3 decreasing for halos as the fraction of 
their mass in substructure increases — since, as we shall 
see, the fraction of halo mass in substructure correlates 
positively with M V i r /M* and consequently inversely with 
halo age. 

We will now address the objects that we intend to re- 
move from the following analysis. These objects have 
abnormally high kinetic energy given their mass. They 
also tend to be lower mass cluster halos (< 10 14 }i~ 1 Mq). 
In Figure IH1 we plot the particle root-mean-square veloc- 
ity dispersion against halo mass. We have ringed all the 
halos that fail the virialization cut. As expected, there 
is clearly a strong correlation between v rms and mass for 
the halos we identify as virialized. The solid line rep- 
resents the relationship between rms velocity dispersion 
and mass for an NFW profile with concentration 6.5 and 
a e( r )/ a i( r ) = 0.5, where a$(r) and cr r (r) are the az- 
imuthal and radial velocity dispe rsions respectively (e.g. 
Eq. 15 of lLokas fc Mamonl {2001)). We choose the above 
concentration as it is approximately the characteristic 
concentration that we obtain for our halo sample (see 
Sec. 15.3(1 . Halos that failed the virial criterion clearly lie 
above the mass - velocity dispersion relation. As they are 
still in the process of relaxing into dynamical equilibrium, 
possibly having recently undergone a major merger, it is 
not surprising that they have an excess of kinetic energy. 



b) 




c) 




Fig. 3. — Projected density of the three example halos (labeled 
a (top), b (middle) and c (bottom)) in Figure 
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Fig. 4. — Scatter plot of (3 vs the virial mass of each halo (includ- 
ing the pressure term). The thick line denotes the median value 
of f3 in logarithmic mass bins. Halos with f3 < —0.2 (circled) are 
not used when calculating the median values in each bin. We take 
M, = 8.0 x lO 14 ^- 1 ^. 
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Fig. 5. — Same as Figure HI but for / s (the fraction of mass in 
substructure) rather than halo mass. The trend-line is calculated 
having omitted halos with (3 < —0.2 (circled). 




Fig. 6. — Scatter plot of halo mass, M vlY against the rms particle 
velocity in each halo. We have circled halos which fail the virial- 
ization cut (i.e. /3 < —0.2). These halos depart from the otherwise 
strong correlation between v rma and M v i r , appearing to have a 
particularly high kinetic energy per particle. The solid line is the 
NFW relation for mass and v rms , assuming a halo concentration 
of 6.5. 



Now that we have defined the halos included in our 
samples, we can proceed to discuss the physical prop- 
erties that we measure. We break these down into two 
subsets: 

• Integral Properties. These are properties that are 
obtained by summing up properties of the individ- 
ual particles in each halo; halo mass, maximum 
circular velocity (as a measure of halo potential) 
and angular momentum. 

• Internal Structure. These describe the distribution 
of the particles within each halo and include con- 
centration, morphology and the distribution of sub- 
structures. 

3.1. Integral Properties 

3.1.1. Maximum Circular Velocity 

In order to measure the maximum circular velocity 
of each halo, we in tegrate the NFW density profile 



ijNavarro et alJ l|1996|) to obtain 



V 2 = tnGrfpi 



In 



Ti + r 



(7) 



where rj and pi = 4p(ri) are the scale radius and density, 
and fit to the circular vel ocity profile of each halo. This 
procedure is exploited bv lKlvpin et al.1 1^)99) to remove 
energetically unbound particles by calculating the escape 
velocity at the radius of each particle relative to the halo 
centre of mass. 

The circular velocity profiles for each halo in our halo 
samples were measured by radially binning particles into 
bins of equal particle numbers. This not only ensures 
that the Poisson errors associated with each bin are 
equal, but also that the radial density of bins is propor- 
tional to the actual density of the halo at any particular 
radius. Using the total number of particles within the 
radius of the outermost particle in each bin, the circular 
velocity at that radius can then be calculated. 

We find that the fitted parameters do not vary sig- 
nificantly with the number of bins used. Furthermore, 
binning weakens the effect of substructure on the profile 
parameters. To quantify the goodness of fit of the NFW 
for each halo, we adopt the relative standard deviation 



0~rel 



'E-=i[(K(r 



V c , n Urj))/V c {r 3 )f 



'3' 



(8) 



and 



where = GM{< Tj)/rj for bin radius 
K.nfwOj) is the predicted NFW circular velocity at this 
radius. The r e has been much recent d iscussion (e.g. 
iKlvpin et aT] l)2001|) . IPower et alJ l)2003|) . lHavashi et alJ 
regarding the conditions that must be satisfied 
in order to ensure numerical convergence at the inner- 
most radius to whic h the circular veloci ty profiles are 
computed. Following Havas hi et alJ (|2004fl , we adopt the 
following three criteria: 

• The number of particles enclosed within the inner- 
most radius, r- m , is sufficient to ensure that five 
times the local collisional relaxation timescale is 
longer than the age of the universe, which is of the 
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order of the orbi tal timescale at the viria l radius 
(see for example IB innev fc Tremainelll987l) 

treiax(r in ) _ N(r in ) r in /V c (r in ) 
) 8ln(N(r in )) r 200 /V 200 




N{r in ) 



ln(N{r in )) 



Pcrit 



-1/2 



> 5 



0) 



where N(ri n ) is the number of particles and p(ri n ) 
the mean density within the radius r*j n . It takes at 
least 15 relaxation times for core collapse to occur 
and 200 for the system to evaporate, hence we have 
chosen 5 as a convergence criterion to allow us to 
probe the inner regions of each halo, without being 
affected by numerical instabilities. 

• The innermost radius is greater than 4 softening 
lengths (e = 3.2h kpc in this simulation) 

• There at least 100 particles within the innermost 
radius 



As in Havas hTet alJ (|2004Tl . we find that the first con- 
straint is the strictest for the vast majority of our halos, 
and virtually every halo that satisfies the third constraint 
automatically satisfies the second as well. 

3.1.2. Spin 

The angular momentum of simulated dark matter ha- 
los is typically been parameterized b y the dimensionless 
spin parameter l|Peebleslll96llT98rl 



A 



_jVe_ 

GM 5 / 2 



(10) 



where M and E are the total mass and energy of the sys- 
tem, G is the gravitational constant and J is the total 
angular momentum of the halo. This definition of spin 
is directly related to the ratio of the centrifugal to grav- 
itational acceleration at the equator of a homogeneous 
rotating body, or the amount of angular momentum re- 
q uired for rotati onal support. 

iPeeblesI (jl969j) showed that the amount by which a 
uniformly rotating_polytrope is deformed, or how oblate 
it is, is related to A in the following way: 



J 2 \E\ 
G 2 M 5 



A 



(11) 



i.e. the flattening of a rotating sphere is directly pro- 
portional to the square of its associated spin parameter, 
where a is the major and c the minor principle axis of 
the halo. 

As we outlined in Section 1, the total energy compo- 
nent of A is difficult to compute due to the necessity 
of accounting for the surface pressure at the outer ra- 
dius of each halo. This is especially the case for ha- 
los in high density re gions of the simulation. Therefore, 
iBullock et alJ (|2001aTl introduced an alternative param- 
eterization, 

A = ^ , (12) 

V2MVR 

which sidesteps this difficulty by removing the energy 
dependence, thus making it more practical to calculate 
from simulated data. We adopt this version of the spin 



parameter, using the virial mass and radius for M and 
R respectively and calculating the circular velocity com- 
ponent, V 2 — GM(R VU ) / R vlr . This retains some con- 
sistency between the two approaches, as A reduces to 
A when measured at the virial radius of a trun cated 
isothermal sphere. IHetznecker fc Burkerti ((2005) have 
performed a comprehensive study of the differences be- 
tween the two definitions of halo spin. In agreement with 
IBullock et alJ l|2001aj) they find there to be little differ- 
ence in the overall distributions of both A and A at z=0. 
However, they also demonstrate that whilst A (A' in their 
notation) does not vary noticeably with redshift, A in- 
creases significantly as redshift decreases. This discrep- 
ancy appears to be due to the impact of the continuous 
acquisition of mass through accretion, which increases A 
but slowly decreases A. 

3.2. Internal Structure Properties 

3.2.1. Substructure 

The procedure we employ to identify substructure and 
build the associated family trees for each cluster halo is 
briefly described in Sect ion 121 The full procedure is com- 
pletely described in lWeller et alJ l|2005j) . One important 
feature of this algorithm is that near neighbor subhalos, 
which are gravitationally bound to one another, are asso- 
ciated together. If this is not done, the subhalo distribu- 
tion might depend on the geometrical smoothing lengths 
used in the definitions of subhalos. Here we will present 
the subhalo mass-function for the 2159 cluster halos in 
our sample, including those tagged as 'unviriliased' . 

3.2.2. Concentration 

Halo concentration is normally defined as the ratio of 
some measure of the outer radius of a halo to an inner 
'scale' radius, normally chosen to be that of the NFW 
profile, i.e. when the logarithmic slope of the density 
profile = -2. Previous studies have used various defini- 
tions of the outer radius of a halo: normally it is chosen 
to be the radius at which the halo density becomes some 
factor of the critical or mean background density of the 
universe at the appropriate redshift. In this study we 
use the virial radius to determine halo concentration, so 
c = Ryir/ri. When comparing our results with studies 
that have adopted an alternative definition of halo con- 
centration (e.g. relative to the mean background, rather 
than the critical density) we use a conversion factor to 
correct for the difference. 

3.2.3. Morphology 

In order to measure halo morphology, we calculate the 
moment of inertia tensor for each halo, Iij = m p J2k r k r ki 
summed over all particles in the halo, where m p is the 
particle mass in the simulation and r\ is the x, y or z 
component of the distance of the k-th particles from the 
cluster center of mass. We then find the magnitude and 
orientation of the principle axes for each halo by diag- 
onalising 1^ and calculating the associated eigenvalues 
and eigenvectors respecti vely. We define the axis ratios 
as in lWarren et al.l l)1992f) . denoting a as the longest axis, 
b the intermediate axis, and c the shortest axis of the mo- 
ment of inertia ellipsoid. 

Previous studies have calculated the inertia tensor 
using an iterative scheme in which only particles in 
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Fig. 7. — Values of Ao and <r as in the literature (see text). 
The large black solid circle is the result of the analysis presented 
here for our main halo sample, the smaller black solid circle for our 
'unvirialised' sample. 



a particular sphere or shell are used to calculate Jy. 
This is then repeated, using the ellipsoid or shell deter- 
mined by the principle axis vectors from the diagonal- 
ized inertia tensor, until until the proced ure converges 
llWarren et alJ Il992|). In agre ement with I.Ting fc Sutol 
( 2002) and iBailin k. Steinmetd l)2004|) we find that this 
approach often does not converge in high resolution simu- 
lations where the over-merging problem — which tends to 
erase subs tructure in low res olution simulations — does 
not apply. Uing &: Sutol l|2002j) adopt a novel approach in 
which they fit halos shapes to isodensity contours. How- 
ever, lAllgood et alJ l)2005[) find that this has the effect 
of suppressing the impact of large substructures near the 
halo core, making higher mass halos appear more prolate 
than they really are. 

4. DISTRIBUTION OF HALO PROPERTIES 

In this Section we describe the distributions of the halo 
properties listed above for the halos in our sample. Un- 
less specifically stated otherwise, we plot separately the 
properties of the halos tagged as 'virialized' and, where 
relevant, those tagged as 'unvirialized'. 



As has been 



4.1. Spin 
established 



in earlier stud- 



ies ( Barnes fc Efstathioil 119871 iWarren et al l 
19921 Ivan den Boschl 119981: Bullock et all 12001a!: 



van den Bosch et al. 2002; Gardner 2 0011 IVitvitska et al.l 



20021 IPeirani et alJ 12004^ . the distribution of A is well 
fit by a log-normal distribution 



P(A) = 



1 



AV27TCT 



exp 



ln 2 (A/A ) 
2a 2 



(13) 



In their analysis, iBullock et alJ l)2001aj) find (Ao, a) 
= (0.042 , 0.5) ( blue, solid square in Figure 0), whilst 
IGardneJ (|2001f) find (0.044,0.5) (red, upward-pointing 
triangle) for the halos that have undergone major merg- 
ers during their formation, and (0.035,0.596) (red cross) 
for th ose that have g rown through accretion. In a similar 
study. IPeirani et alJ l|200^) found (0.036,0.57) (blue, five- 
pointed star) over all the halos in their sample, whilst Ao 
ss 0.044 (filled blue diamond) for the merger halos and Ao 
pa 0.031 (blue diamond) for the a ccretion dominated ha- 
los (see their Figure 7) . Likewise, iHetznecker fc BurkertJ 



(2005) obtain (0.039,0.53) for their entire halo sam- 
ple (black plus sign), (0.041,0.53) for their merger ha- 
los (black, downward pointing triangle) and (0.037,0.55) 
for the accretion dominate d halos (filled, b lack d own- 
ward pointing triangle). iVitvitska et al.l l)2002|) find 
(0.0445,0.563) from their numerical simulations (green 
circle) and (0.0465,0.556) (green square) from their ran- 
dom walk model (based on the simulations). Whilst 
stud ying the dependence of A on environmental den- 
sity, lAvila-Reese et all (|2005|) obtain (0.032, 0.59) for 
their 'Void' halo sample (open, red six-pointed star) and 
(0.034,0.60) for their 'Field' halo sample (filled, red six- 
pointed star). Finally, by modeling the acquisition of 
spin by halos from both the orbital angular momentum 
of merging satellites and through the effects of a large- 
scale tidal fi eld (tid al tor que theory — see for example 
!Peeblesil969l lWhitdll98l . iMaller et"afl pOOl find that 
both scenarios produce a log-normal distribution of halo 
spin with « (0.04,0.5) (filled magenta triangle). Overall, 
most recent studies fix Ao in the range 0.031-0.045, with 
dispersions varying between 0.5-0.6. 

The value we obtain for Ao (0.033 — see lower panel 
in Figure |H1 and the large filled black solid circle in Fig- 
ure Q) is in broad agreement with these studies, how- 
ever we obtain a slightly broader distribution with a 
= 0.645. Fo r comparison , we al so plot the distribution 
obtained by IMaller et all (|2002) (dotted line in Figure 
|SJ). Interestingly, our results appear t o be very simi- 
lar to those obtained bvlGardnerl (|2001|) , IVitvitska et alJ 
( 2002) and IPeirani et alJ (|2004l) for the halos in their 
sample that have grown through accretion only, having 
not undergone any ma jor mergers during th eir formation. 
IVitvitska etldl (|2002|) and lGardnerl l)2001|) both suggest 
that increased merging activity will result in a higher 
halo angular mome ntum — a result that is verified by 
IPeirani et al. ( 2004) , who monitor angular momentum 
as halos form in their simulations. However, the higher 
dispersion that we measure than any previous study does 
indicate that there are still a significant number of high 
Ao halos in our sample, indeed, over 300 halos have Ao > 
0.06. 

The top panel in Figure OH describes the spin distri- 
bution for halos that we have tagged as 'unvirialized', 
using the criterion outlined in Sec. 12.21 This is also well 
described by Eq. El with A = 0.049 and a = 0.612 
(small, black filled circle in Figure 0). A larger value of 
Ao is consistent with the idea that these halos have not 
yet reached a state of dynamical equilibrium and have 
greater kinetic energy (partly in the form of angular mo- 
mentum) than their virialized counterparts. 

4.2. Mass function and fraction of mass in substructure 

The bottom panel in Figure El displays the distribution 
of the fraction of mass contained in substructure for our 
main halo sample. The median value of f s is 0.056, where 
f s ranges between < f s < 0.65. Note also that there 
is a large tail to high values of f s . The upper panel 
of Figure |S] demonstrates that the unvirialised halos do 
have on average more substructure, with a median value 
of f s = 0.196. This indicates that these are likely to 
be young clusters in which there has not been sufficient 
time for the halo core to remove a significa nt fraction 
of subh alo mass through tidal ab lation. Both lGill et alJ 
()2004cj) and lDe Lucia et alJ l)2004|) , who resimulate 8 and 
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Fig. 8. — (bottom panel) Di strib ution of the halo spin parameter, 
A (crosses) and fit of Eg uation ll3l f solid line), where Ao = 0.033 and 
o~ = 0.645. F or comparison, we plot a distribution with Ao = 0.04 
and <t = 0.5 I Mailer et al. 2002). (top panel) Same, but for those 
tagged as 'unvirialized', where Ao = 0.049 and a = 0.612. 



11 cluster halos respectively, and lGao et alJ ((2004), who 
resimulate halos over a wider range of masses, obtain 
mean values for f s close to the ones presented here. 

In the distribution in the lower panel of Figure ED there 
is tail of high f s halos. These are halos that have very 
recently undergone a major merger (although they have 
passed our virialisation criterion) and have a dual core 
at their center — the mother halo and a slightly less mas- 
sive subhalo. In the upper panel of Figure ED we plot the 
mass of the largest subhalo as a fraction of the total mass 
in substructure against f s . It is clear that, on average, 
halos with a large fraction of mass in substructure host 
a subhalo of almost equivalent mass to the mother. This 
is also demonstrated in the lower panel, where we plot 
the distance between the center of mass of the halo and 
its position of peak density, as a fraction of the virial ra- 
dius, against f s . The density peak, by definition, resides 
in the mother halo. However, for halos with an especially 
massive subhalo, the center of mass will lie somewhere 
in between the density peak of the mother and that of 
the subhalo. For these halos, an NFW profile provides a 
poor fit (see Figure IT5|l to the radial density profile due 
to the impact of the large subhalo. The correlation indi- 
cates that the center of mass moves further from the peak 
density position as f s increases — suggesting that much 
of the substructure is clustered on one side of the halo. 
This can be explained if halos with large f s are more 
likely to have one large piece of substructure, rather than 
lots of smaller subhalos that just happen to be located 
on one side of the mother halo. 

In Figure ITT1 we plot the mass distribution of the sub- 
halos we obtain from our sample of 2159 clusters, includ- 
ing those tagged as 'unvirialised'. The dotted histogram 
is the distribution for all the sample halos, while the 
solid histogram only includes the 2085 virialized halos. 
The number density iV su b of subhalos of mass Af sll h can 
be described by a Sch echter function ijPress fc Schechterl 
Il97i ISchecMerllT97^ : 



dK 



sub 



= Nj -^L cxp -^0. (14) 
dln(M sub ) \M* sub J [ M*sub_ 

which is a power law with an exponential cut-off at large 
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Fig. 9. — (bottom) Distribution of the fraction of mass in sub- 
structure for our main halo sample, (top) The equivalent distribu- 
tion for those halos that were tagged as 'unvirialised' 
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Fig. 10. — (top) Ratio of the mass in the most massive subhalo 
to the total mass in substructure, as a function of the total fraction 
of mass in substructure, / s . The plot demonstrates that for halos 
with f s > 0.3 the substructure is largely dominated by a single 
massive subhalo of mass comparable to that of the mother halo. 
(bottom) Scatter plot of the distance between the point of peak 
density and the center of mass, as a fraction of the virial radius, 
for each halo against f s . 



mass scales, where N* is a normalization constant. In 
Figure ITT1 we show the mass function of subhalos (solid 
line) where the substructure mass is expressed in units 
of the halo mass; this brings potentially different mass 
halos onto an equal footing. At a mass scale of about 
Msub = 3 x 10~ 3 Mh the function drops, because at this 
stage our finite sample size and the lower mass limit of 
our identified halos of 7.6 x 10 10 h^ 1 MQ truncate the dis- 
tribution. 

We fit this distribution up to the truncation point with 
the Schechter function (Equation I14fl . varying all three 
parameters, and obtain a — 0.5 and M„ su (, = 0.3-M/j 
(solid line). However, there are large correlations or pa- 
rameter degeneracies amongst a, N*, and M* su b, which 
lead to marginalized errorbars in a of ±0.2 and in M* su }, 
of ±0.7Mft. In order to compare with previous work we 
omit the cut-off, M* su (,, and just fit for the slope and am- 
plitude. In this case we obtain a = 0.75 (dashed line). 
Still there is a remaining correlation between a and 
and the marginalized errorbar on the slope is ±0.12. If 
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Fig. 11. — Mass distribution of subhalos from the 2159 most 
massive halos in our simulation as a function of the ratio of subhalo 
mass to cluster mass. The dotted histogram is the distribution of 
the complete sample, the solid histogram for the 2085 virialised 
halos only. The solid line is the Schechter function fits as discussed 
in the text. The dashed line is obtained for just fitting the slope. 
In both distributions halos with less then 30 particles have been 
dissolved into their mothers, thus setting a lower mass limit for the 
subhalos in our sample. 



we just fit for the slope we obtain a — 0.75 ± 0.12. 

Furthermore, we find the distribut ion becomes flat be- 
low M sllh » 3 x lQ- 3 M h . Others iGhigna et all 1200(1 
Springel et all 120011 IDe Lucia et all l2004t lOao et all 
2004) have found similar results with regard to the power 
law component; but we have added an exponential cut- 
off at the high mass end, which is only possible of our 
large halo sample. We note the cutoff in high mass ob- 
jects away from the simple pow er law, as is pr edicted 
from the Schechter prescription (|Schechterlll9 76') . Such 
a cutoff is of course essential, as there cannot, by defini- 
tion, exist subhalos more massive than the mother halo. 
At this stage the distribution of subhalos is steeper than 
the global distribut ion of all the halos in our simulation 
fWeller et al.l l2005") by almost 0.2 in the exponent. Our 
explanation for this is that while the original input distri- 
bution to a growing halo must correspond to the global 
average distribution, the more massive subhalos are more 
quickly destroyed as dynamical friction subjects them to 
more rapid tidal ablation, thus steepening the profile. 

Note that the slopes of the Schechter fu nction obtained 
here a re lower than the ones obtained in IDe Lucia et all 
|200l . which are in the range a ~ 0.85 — 1.13. How- 
ever if we fit for the slope alone the 2-cr region of 
our fit i s 0.51 < a < 0.99, which overlaps with this 
range of IDe Lucia et all l)2004lL T he difference between 
our result and previous discussions iGhigna et al.lll 998b; 
Moore et al1ll999t ROvpin et al.lll999l: iHelmi et al.ll2002 : 
De Lucia et alJl2004l) is due to the fact that we include 
all halos, regardless of their morphology, while previ- 
ous studies picked halos which are not in merger states. 
While the choice of relaxed halos clearly produces a bias, 
our sample might be biased as well, since one expects 
more mergers among the most massive halos in the sim- 
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Fig. 12. — (bottom) Distribution of halo concentrations, includ- 
ing (crosses) and omitting (circles) halos with f s > 0.3. Solid and 
dashed lines are the best fit log-normal distributions, (top) The 
distribution for halos ta gged as 'unvirialised' according to our viri- 
alisation criterion (Sec. 12:21 . 



ulation than in a randomly selected sample. 

4.3. Concentration 

It has be e n dem onstrated by I.Iind (|2Q0Q) and 
iBullock et alJ (|2001bf l that the distribution of halo con- 
centrations at any particular redshift can be well de- 
scribe d by a log-n o rmal distribution, of dispersion a « 
0.22. IDolae et all l|2004() showed that this result is in- 
dependent of the adopted cosmological model. It is also 
thought that the dispersion in concentrati ons is related to 
the d ispersion of halo formation epochs (Wcchslcr et al. 
2002). 

Figure ^| shows the best fit log-normal distribution 
(solid line) to our halo sample. At first glance, this Figure 
appears to suggest that our halo sample (crosses) does 
not conform to such a distribution. This is mainly due to 
the excess of halos with a very low concentration (< 3). 
These halos appear to have significantly more mass in 
their outer regions than is suggested by an NFW profile. 
Indeed, as demonstrated in the upper panel of Figurc lTSl 
the NFW profile provides a very poor fit to these halos; 
the reason for this is that these halos have a considerable 
fraction of their mass in substructure (Figure 1131 lower 
panel). Halos with a concentration of less than ^2.1 are 
those that have a massive (bound) daughter at or near its 
virial radius. For these halos, the circular velocity profile 
continues to rise beyond i? V j r . The open circles in Figure 
1 121 represent the distribution of halo concentrations once 
all halos with f s > 0.3 are removed. It is clear that this 
sub-sample is well described by a log-normal distribution 
(dashed line) with a dispersion a c = 0.22. 

The top panel displays the distribution of concentra- 
tions for the halos tagged as 'unvirialized'. The distri- 
bution is also well described by a lognormal distribution, 
but peaks at a far lower concentration than the main 
halo sample. In these halos, much of the matter that 
will eventually be accreted onto the halo core is still in 
the outer regions, thus reducing the concentration. Fur- 
thermore, they also have (on average) a higher fraction 
of their mass contained within substructures, many of 
which are located outside the central regions of the halo. 
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Fig. 13. — Scatter plots of halo concentration (upper panel) and 
the fraction of mass in substructure (lower) with respect to the 
goodness of the NFW profile fit (as defined by Equation 151. 
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Fig. 14. — (Bottom:) histogram of the distributions of the halo 
axial ratios, b/a (solid line), c/b (dashed) and c/a (dotted). The 
median values agree quite well with previous studies. ( Top:) same, 
but for the halos tagged as 'unvirialised'. 



4.4. Morphology 

Figure n~T1 displays the distributions for each principle 
axis ratio. The ratio of the intermediate to minor axes 
is in general greater than that of the major to interme- 
diate axes. This implies that the halos are more prolate 
than oblate — a resu lt found by se veral previous stud- 
ies (IDavis et al.ll 198511 Warren et ai]|1992 ICole fc Lacevl 
119961 IKasun fe Evrardl 120051: IBailin fc Steinmetd I2004J) . 
None of the three distributions resemble a normal distri- 
bution, unlike those found bv Uing fc Sutol J2002), wn0 
resimulate a small number of halos over a broad mass 
range (10 1 2 - 10 14 M m ). A s both IBailin fc Steinmetd 
(2004J and|K asun fc Evrardl l)2005|) have found, we see 
a tail towards very low axis ratios, corresponding to ex- 
tremely flattened halos. 

The mean values we obtain for the axis ratios are b/a 
= 0.817 ± 0.098, c/b=0.867 ± 0.067 and c/a = 0.707 
± 0.095 respectively. Again, the se values agree quit e 
well with the axis r atios found by |Warren et_aL| l)1992|) ; 
ICole fc Lacevl 1I1996D andlKasun fc Evrardl <|2005^ . In the 
case of IWarren et al J l|1992j) we find our results match 
those they obtain for the largest ra dius used in their iter- 
ative ellipsoidal fitting technique. iDubinski fc Carlberd 



l)1991[) obtain a very low mean c/a ratio of 0.5 for the 
isolat ed collapse of halos wi th lOOkpc rad ius. However 
both IDubinski fc Caxlberel lj!991|) and IWarren et all 
( 1992) also find that the halos become steadily more 
spherical as the length of the semi-major axis of the shell 
was increased in their analy sis. From th e Gaus sian fit 
to their distribution of c/a, Uing fc Sutol (|2002f) find a 
mean value of 0.54 with a dispersion of 0.113. They then 
construct a conditional probability distribution for b/a 
over certain ranges of c/a, yielding a maximum likeli- 
hood value of 0.77 for b/a , given their mean value of c/a. 
IBailin fc SteinmeO l)2004|) break each of their halos down 
into concentric spheres before measuring the inertia ten- 
sor in each shell (rather than for the entire halo) so it is 
difficult to compare directly the values we obtain for the 
mean axis ratios. Comparing to their most massive halos 
(10 13 — 3 x 10 14 M©) we find that the axis ratios we obtain 
are greater (i.e. closer to 1), with the exception of those 
they obtain using only the inner most shells (se e their 
Figure 3). Finally, recent results bv lAllgood et al.1 ( 2005) 
measure a c/a ratios of 0.6 ± 0.1 for their galaxy mass 
halos. These results are a lso in reasonabl e agre ement 
with observational work by Basilako s et alJ ((2000), who 
measure the projected cluster ellipticities of APM galaxy 
clusters based on moments of the smoothed galaxy dis- 
tribution. They find that a prolate spheroidal model fits 
the data best and (without excluding clusters that show 
significant substructure) measure a c/a axis ratio of 0.65 
± 0.15. 

Overall, the halos in our simulations are slightly more 
spherical than those in previous studies. This is possibly 
due to our selection policy of removing from our analysis 
subhalos that have their center of mass just outside of 
the virial radius of the halo. Large pieces of substructure 
located just outside of the halo would, if included as part 
of the halo, have the effect of making it appear more 
prolate and thus decreasing the c/a axis ratio. In order to 
test this we have repeated the analysis, but now including 
all substructure that has at least one particle inside of 
the virial radius of the halo. Figure ITol compares the c/a 
distributions obtained from this new dataset with the 
original. It is clear that the new histogram has many 
more halos with a low c/a axis ratio than the original, 
especially in the tail. The new mean values for the axis 
ratios are b/a = 0.79 ± 0.121, c/b = 0.859 ± 0.071 and 
c/a = 0.678 ± 0.113. The fact that the mean minor-to- 
intermediate axis ratio has decreased the least suggests 
that the halos have become more prolate with the added 
substructure. 

The top panel in Figure ^] shows the distribution of 
the axis ratios for the halos that we have tagged as being 
'unvirialized' according to our criterion outlined in Sec. 
12.21 They indicate that these halos arc less spherical 
than their more relaxed counterparts, probably due to 
their larger fraction of mass contained in substructure. 
We investigate further the impact of substructure on halo 
morphology in Section H3 

Another way in which the shape of simulated dark mat- 
ter ]mk2s_£anj3eck^^ 

T l|Warren et al.ll992tlBlm"in fc Steinmetzl2004|) . defined 
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Fig. 15. — Solid histogram: distribution of the c/a axis ratio. 
Dotted histogram: the same, but substructure just outside of the 
virial radius is included in the calculation of the inertia tensor for 
each halo. The solid and dotted arrows indicate the median values 
of each distribution. There is a clear shift towards lower values of 
c/a when including outlying substructure, suggesting a slightly less 
spherical morphology. 
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Fig. 17. — Maximum circular velocity as a function of halo 
mass, having removed halos for which the NFW profile provides 
a poor fit due to large amounts of substructure. Straight line is 
best fit of function Vcmax = A(M vil /M*) a , with a = 0.31 ± 0.08. 
The dashed line is the NFW prediction based on our me asur ed 
correlation between halo mass and concentration (Equation 1161 . 



30 
20 
10 


300 
250 
200 
150 
100 
50 
0, 



P < -0.2 



P>-0.2 
Oblate 



Triaxial 



Prolate 



0.2 



0.4 0.6 
Triaxiality 



Fig. 16. — Distribution of the triaxiality, T (see Equation ll5l . for 
our main halo sample (bottom), and those tagged as 'unvirialised' 
(top). Thick black lines indicate the medians of each distribution. 
Both panels show that there is a clear preference for halos to adopt 
a more prolate, rather than oblate, morphology. 



The full range of triaxiality falls into three regimes: 
oblate (or 'pancake-like') halos correspond to T < 1/3, 
prolate (or 'cigar- like) to T > 2/3 and fully triaxial ha- 
los to the intermediate range. From Figure ^3 it is clear 
that the majority (52%) of halos lie in the T>2/3 range 
and are therefore more prolate than oblate, with the dis- 
tribution peaking at T~0.75. A significant number of 
halos (38%) are completely triaxial — the principle axis 
lengths are distributed over a significant range. 

5. CORRELATIONS WITH MASS 

In this Section we investigate how the halo properties 
we measure correlate with halo mass. This is important 
for determining whether certain halo characteristics are 
self-similar, or whether they vary significantly depending 
on the age or phase of evolution of a halo. All the results 
presented in this Section are for the 'virialised' halos in 
our sample only. As most previous studies have found 
mavis et a,l.ll1985UWarren et a,l.ll1 992HAvila-Reese et a.1.1 
2005), we see no measurable correlation between halo 



tion 



mass and spin. 

5.1. Maximum Circular Velocity 

Figure 1171 plots maximum circular velocity as a func- 
tion of virial mass. We have removed all halos with 
f s > 0.3 as the NFW profile provides a poor fit to 
many of these halos due to their large quantity of sub- 
structure (see Figure I13J) . Fitting a simple power-law 
Vcmax oc M 1 /", we measure a = 0.31 ± 0.08 (solid 
lin e). This is very close t o the relation s hip obt ained 
bv iKravtsov et alJ l|2004rfl . iBullock et all l|2001bD . and 
H^v^gh^g^alJ 12QQ3J) (for their subhalos in the latter). 
Bulloc k et alJ l)2001b|) explain the departure from that 
predicted by standard scaling of the virial parameters, 
M vir cx V^ r , through the observed mass dependency 
of the concentration parameter — a result that we also 
find (Section OJ). For the NFW profile, V cmax cx 

^ 3 /(c), where /(c) is a function of halo concentra- 
which in turn de pends on M„ r (see, for exam- 
ple, lAvila-Reese et all (pH). In Section E3 we find 
a power-law relation between c and M V i r . The dashed 
line in Figure 1171 represents the corresponding relation- 
s hip between M„,> and V rrn ^ ■ 

IKravtsov et alJ l)2004b|) follow the evolution of dwarf 
halos as they first accrete mass in isolation and then un- 
dergo significant mass loss through tidal stripping whilst 
being captured by more massive neighboring halos. They 
find the M-V^ max relationship holds throughout the dif- 
fering phases of subhalo evolution. Hence, it appears that 
this correlation is independent of halo mass and forma- 
tion epoch. This helps to explain the tightness of the 
relation in Figure IT71 

The maximum circular velocity of a halo can be in- 
terpreted as a tracer for the depth of halo central po- 
tential. It is more reliable for this purpose than halo 
mass, which can be difficult to define, and especially so 
for subhalos. Figure ITS1 compares V c 2 max and M v i r as a 
function of 4> c , the gravitational potential at the point of 
peak density in the halo (as identified from the smoothed 
particl e distribution by th e Denmax routine in our halo 
finder ijWeller et alJl2005j) 1. The solid lines in each plot 
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Fig. 18. — Maximum circular velocity (top) and halo mass (bot- 
tom) as a function of halo central potential, <p c . Lines are best fit 
power laws X oc rf>2, where X is V^ max and 7 = 1 in the top panel 
and M v ; r and 7 = 1.5 in the bottom. 
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Fig. 19. — The fraction of mass in substructure as a function of 
halo mass. Points with errorbars represent the median values of f B 
in logarithmically spaced mass bins. The dashed line is the best fit 
power law, with slope 0.44±0.06. We take M* = 8.0x 10 14 /i -1 A/q . 



represent the best fit power-law. We obtain V^ max oc 4> in 
the upper panel and M vir /M r oc </> 15 for the lower. Al- 
though both plots show clear correlations, that of V cmax 
and 4> c (top panel) is clearly the tighter. To quantify 
this, we adapt Equation|Hlto measure the goodness of fit 
of the power-law in each case, obtaining <j re i — 0.08 for 
the Vcmax ~ 4>- correlation and a re i — 0.19 for the mass— cA 
correlation. This confirms that V cma x is the more robust 
measure of halo central potential. 

5.2. Substructure 

We plot in Figure ITT)1 the fraction of mass in substruc- 
ture as a function of halo mass. The dashed line repre- 
sents the best fit power-law, f s — A(M V i r / 'Af*) m , with 
slope 0.44 ± 0.06 and amplitude 0.14 ± 0.02. Hence, 
higher mass halos tend to have contain a higher frac- 
tion of their mass in susbtructure. This is an im- 
portant result, demonstrating that cluster-mass halos 
are not self-similar. Several studies of small sam- 
ples of halos selected from lower resolution simulations 
have suggested the dark matter halos are self-similar 
in terms of their subhalo populations: low mass ha- 
los appear lik e 'rescaled versions' of hi gher mass halos 
(iMoore et alJ 119991: iGhigna et all l200Gt iDe Lucia et all 



Jm 

he 



120041 iReed et al.ll2005blh With a statistically signifi- 
cant sample of well reso l ved h alos, we find this not to 
be the case. iGao et all ()2004|) . who compare the sub- 
halo populations of a small sample of halos over a wide 
mass r ange, found a similar result to the one presented 
here. Ivan den Bosch et al.l l)2005j) constructed a semi- 
analytical model to calculate the subhalo mass func- 
tion, adjus ting the fr e e par ameters to match the sim- 
ulations of IGao et all l)2004|) . iTormen et all <|2004|) and 
iDe Lucia et alJ |2004|) . They found that the slope and 
normalization of the mass function to be dependent on 
the ratio of halo mass to t he ch aracte ristic non-linear 
mass scale. iZentner et alJ l|2005|) and iTavlor fc~ Babul 
( 2004) have also proposed models with similar proper- 



ties, arguing that as higher mass halos have formed more 
recently, they should contain a higher fraction of their 
mass in substructure. 

5.3. Concentration 



In their early studies, NFW postulated that the con- 
centration of a halo is an indicator of the mean density 
of the universe at the time of its collapse ijNavarro et alJ 
1997). Hence, to conform with the current paradigm 
of hierarchical structure formation, halo concentra- 
tion ought to d ecrease with increas i ng mass. Sub- 
sequent studies JBullock et al.l l2001bt lEke et all 120011: 
Dola g~ alJl2004l) have suggested analytic formulas pre- 
dictin g halo concentratio n as a function of mass and red- 
shift. iDolag et alJ {2004) suggest a simple fitting formula 



c(M, z) 



1 + z 



M 



(16) 



finding a = -0.102 and c = 5.97 for their ACDM 
simulation, with a almost independent of cosmology. 
The cor responding parameters obtained by Bullo ck et alJ 
<|2001bJ) are a = -0.13 and c = 5.60. (It should be 
noted that both these studies calculate the concentra- 
tion by defining the virial radius as enclosing an over- 
density of A c times the mean background density rather 
than the the critical density as we do in this analysis. We 
therefore have applied the appropriate correction to their 
resu lts to enable comp a risons with our own.) Addition- 
ally, lAvila-Reese et all ()2005|) found a logarithmic slope 
of —0.14 for c(M V i r ), for the halos in high density regions 
in their simulati ons. Analysing halo concentrations at 
sev eral redshifts.lBullock et, aU pOOlhfl , lEke et, all <|200lD 
and lNavarro et al.l l)2004[) all show that the concentration 
parameter for a given mass is lower at higher redshifts. 
By following th e mass accr etion history of simulat ed ha- 
los, [Zh^iLaD (mil) and lTasitsiomi et alJ ((2004) both 
find that halo concentration increases during the slow ac- 
cretion phase of halo formation, having remained fairly 
constant during an initial merger period. 

The median trend (points with error bars) in Figurel2*Ul 
clearly demonstrates the predicted anti-correlation be- 
tween concentration and virial mass at a constant red- 
shift for the halos in our sample. As in Section I5TT1 we 
have removed all halos with f s > 0.3 as the NFW profile 
provides a poor fit to many of these halos due to their 
large quantity of su bstructure. T he solid line represents 
the best fit for the IDolag et alJ l)2004|) fitting formula, 
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Fig. 20. — Scatter plot of concentration as a function of virial 
mass, having removed halos for which the NFW profile provides 
a poor fit due to large amounts of substructure. The points with 
error bars are the median conce ntra tion in logarithmic mass bins, 
the solid line is best fit of Eq. 1161 to the data (a = —0.12 and 
c = 6.47). 



with a = -0.12 ± 0.03 and c = 6.47 ± 0.03. However, 
there is also a considerable amount of scatter in the dis- 
tribution of points making it difficult to accurately deter- 
mine the slope of the co rrelation. In their recent study, 
lAvila-Reese et all ll2005|) showed that halo concentration 
is also dependent on the local environmental density: ha- 
los in dense regions are on average more concentrated 
and have higher centra l densities. As demonstrated by 
ISheth fc Torrneil (pOO^ . halos in dense regions accumu- 
late their mass earlier than their counterparts in lower 
density regions. Hence, low mass halos forming in low 
density regions will have a lower concentration than pre- 
dicted by the lDolae et alJ J2004J) M vir — c relation. This 
can be seen in Figure [201 — there is a significant tail of 
low concentration halos at lower masses. 

5.4. Morphology 

In order to investigate how halo morphology is de- 
pendent on mass, we plot in the upper of Figure CHI 
the minor to major axis ratio versus halo mass. The 
points with errors bars mark the median trend in log- 
linear mass bins. The solid line indicates the best fit 
power law, c/a oc {M mr /M«) s , with S = -0.049 ± 0.007. 
The slope, S, that we obtain is in very good agreement 
with a comprehensiv e recent study of halo morphology by 
lAllgood et all l|2005|) . who measure S = -0.050 ± 0.003 
over a wide range of h alo masses. In another recent 
paper, iPaz et af] lj2005f) measure a steeper r elationship, 
with d = -0.056 ± 0.003. Ding fc Sutol |2002^ . who anal- 
yse a similar halo mass r ange to ours, find a st eeper slope 
still. As pointed out by Allgoo d et alJ l)2005|) . who per- 
form an in-depth comparison of halo morphologies from 
different simulations and measurin g techniques, this is 
probably due to the approach that Uing k, Sutol (2002) 
adopt of measuring the shapes defined by isodensity con- 
tours, thus ignoring the mass distribution in the central 
regions of their halos. Essentially, this suppresses the ef- 
fect on halo shape of large substructures near the halo 
core. Regardless of the shape measuring technique used, 
it is clear that in general halos become less spherical with 
mass. 
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Fig. 21. — Scatter plot of halo virial mass, M v - lr , versus the minor 
to major axis ratio (top) and halo triaxiality (bottom). The solid 
line in the top panel represents the best fit power law, with slope 
—0.049 ± 0.007 (see text). Points with error bars denote median 
value and standard deviation in log-linear mass bins. The solid 
line in the bottom panel is the median trend line. The plots show 
that higher mass halos are less spherical (decreasing c/a) and more 
prolate (T approaching unity). 



In the lower panel of Figure |^ we plot halo triaxiality 
(as defined in Equation I15JI versus mass. As implied by 
the size of the error bars, there is an enormous amount 
of scatter. However, there is a preference for high mass 
halos to be more prolate. Therefore, the apparent trend 
is for dynamically 'younger' halos (larger M V i r /M*) to be 
more prolate. Low mass halos are older than their higher 
mass counterparts; they have had longer to dynamically 
relax and obtain a spherical morphology than the higher 
mass, more prolate halos. Furthermore, higher mass ha- 
los contain a higher fraction of their mass in substructure, 
which also influences their shape. This is investigated 
further in the following Section. 

6. THE IMPACT OF SUBSTRUCTURE ON HALO 
PROPERTIES 

We now investigate the impact of substructure on the 
physical properties of the dark matter halos. Specifically, 
the purpose is to determine how the distribution of sub- 
halos effect the morphology, concentration and spin of a 
cluster. All the results presented in this Section are for 
the 'virialised' halos in our sample only. 

6.1. Mass distributions and concentration 

In the previous Section we found that halo 
concentration — an indicator of the formation epoch of 
a halo — decreases with halo mass. We find there to be 
an even tighter anti-correlation between halo concentra- 
tion and the total fraction of mass in substructure: halos 
with more substructure have lower concentrations. This 
is evidence of the continuous effect of the tidal strip- 
ping of subhalos. Halos with large amounts of substruc- 
ture formed more recently, when the background den- 
sity of the universe — and therefore central densities and 
concentrations — are lower. The less dense outer parts 
of infalling subhalos are stripped off in the outer regions 
of the halo, and the denser core spirals into the central 
region due to dynamical friction, thus increasing the con- 
centration of the halo once the subhalo is fully digested. 
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Fig. 22. — Scatter plot of the fraction of mass in substructure 
against halo triaxiality. Error bars represent the standard deviation 
in each bin. There appears to be a slight trend for triaxiality 
to increase with f s , although the errors on each point are large. 
This indicates a tendency for halos to become more prolate as the 
fraction of their mass in substructure increases. 



6.2. Contribution of substructure to halo morphology 

We have already shown that substructures lying just 
outside of the virial radius have a measurable effect on 
the overall distribution for each axial ratio. Figure E3 
is a scatter plot of the fraction of mass in substruc- 
ture against triaxiality. Again, there is an appreciable 
amount of scatter, although the apparent trend is for 
substructure-rich halos to be more prolate. To further 
quantify the effect of substructure, we recalculate the 
c/a axis ratio for each halo, but having first removed all 
particles that belong to substructure, leaving particles 
that belong only to the 'mother' or 'host' halo. We also 
repeat this exercise, but removing only the most massive 
subhalo in the cluster and retaining the rest. Figure l2~3l 
shows the change in the minor/major axis ratio when 
these two sets of particles are removed. These distribu- 
tions indicate that the particles belonging to the mother 
halo only are more spherically distributed than the rest 
of the halo mass. Therefore, the effect of substructure 
is to make cluster size halos appear more prolate. The 
dashed histogram demonstrates that much of the impact 
of substructure on the minor/major axis ratio is due to 
the most massive subhalo in each halo. Smaller substruc- 
tures, regardless of their position in the cluster, appear 
to have less influence. 

6.3. Effect of substructure on halo spin 

Figure [21 shows the relation between halo angular mo- 
mentum and the fraction of mass in substructure. The 
median value of A increases from 0.030 for halos with 
f s < 0.2, to 0.051 for those with f s > 0.2. This indi- 
cates that halos with higher f s — those more likely to 
have undergone major mergers — have had their angular 
momentum boosted through the acquisition of smaller 
halos. 

To investigate this further, we now consider the contri- 
bution to the total spin by the mother halo only and by 
the substructure. To calculate these quantities we define 
a basic unit angular momentum per unit mass 
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Fig. 23. — Histograms of the c/a axis ratio having included 
(solid) and omitted (dotted) the substructure whilst determining 
Iij for each halo. Also plotted is the distribution having removed 
only the most massive subhalo (dashed). Arrows denote the me- 
dian values of each distribution. The dotted histogram shows that 
when the substructure is discarded, thus leaving only the mother 
halo, the shape becomes more spherical (c/a closer to unity). It 
is clear that the largest subhalo is significant in determining clus- 
ter morphology: removing these results in an obvious shift in the 
overall c/a distribution. 




Fig. 24. — Spin parameter A as a function of the fraction of mass 
in substructure, f s . The median trend line shows there to be, on 
average, a small increase in halo spin, as the fraction of halo mass 
in substructure increases. 



where i? v ir and M v - lr are the total halo energy and mass 
respectively. We then calculate the angular momentum 
per unit mass directly for particles belonging solely to the 
mother halo (A m ) and for those belonging to substructure 
(A s ) in terms of [j]: 



A, 



(J m /M m ) 



(J s /M s ) 
[J] 



(18) 



[j] = [2GR V1I M V 



(17) 



where J m and J s are calculated using particles belonging 
only to the mother halo or substructure, respectively. 
Note if f s = 0, then X s is also zero and X m returns to 
our original definition of A. 

Figure [25] shows the overall distributions obtained for 
X s and A m compared to the original distribution. We 
have also plotted the best fit log-normal curves (Equation 
113(1 where Ao m = 0.029 and a m = 0.608 for the mother 
halo, and Ao s = 0.144 and a s = 0.863 for the substruc- 
ture. These results clearly show that particles belonging 
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Fig. 25. — Distributions of A (circles), A m (crosses) and A s 
(diamonds), and the respective fits to a log -nor mal distribution 
(solid, dashed, and dotted lines; see Equation 1131 . A m and X a are 
calculated using particles belonging either only to the mother halo 
or substructure, respectively. A is plotted on a log scale due to the 
much greater angular momentum per unit mass (see Equation 1181 
of particles belonging to substructure. The fit parameters are Xo m 
= 0.029, cr m = 0.608 for the mother halos only, and A 0s = 0.144 
and cr s = 0.863 for the substructure only. 
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Fig. 26. — {top) Similar to Figure |2*11 except for A is calculated 
only using particles that belong to substructure in each halo. The 
median trend-line shows a weak d ecre ase of A s with f s — the oppo- 
site to the trend demonstrated in !24l for total halo spin, (bottom) 
Same as above except for A calculated only using particles that 
belong to the mother halo. Median trend-line shows an increase of 
A m with f s . Combined, these plots demonstrate that the increase 
of total halo spin with f a is due only to the increase of A m with 



to substructure have a far greater average angular mo- 
mentum per unit mass than those belonging solely to the 
mother halo. Accreted particles will fall into a halo in a 
generally spherically symmetric manner, and will there- 
fore induce only a smaller change to the total angular 
momentum compared to the capture of a subhalo — an 
injection of mass from one particular direction. Impor- 
tantly, in the next Section we find that subhalos tend to 
be found in a prograde orbit with respect to the over- 
all angular momentum of the mother halo, implying a 
preferential direction of infall. However, there is also a 
selection effect here. Subhalos adopting plunging orbits 
directly towards the cluster center of mass, and there- 
fore possessing a lower orbital angular momentum, will 
be more rapidly destroyed through their exposure to tidal 
forces from the halo core. This introduces a bias favoring 
the survival of subhalos in more spherical orbits. 

Overall, we expect the velocity distribution of captured 
dwarf halos to be fairly broad — those that form near to 
their eventual parents will have a lower velocity relative 
to the mother than those that have formed further away 
and have thus gained more kinetic energy as they fell 
into the potential well of their future host. This effect 
explains why <r s , the dispersion in A s , is so great. 

In order to investigate the effect of substructure on 
the spin of the mother halo, we plot \ s and A TO against 
f s in Figure [2lfl In the equivalent plot for the total halo 
spin, we found there to be a correlation between A and f s 
(Figure [21 ■ We find a similar result for the mother halo 
spin, but a weak anti-correlation between substructure 
spin and f s . At first glance this may seem strange, but 
it is merely due to the fact that halos with a large f s tend 
to contain at least one very massive piece of substructure. 
Being located near the halo center, the subhalo does not 
contribute as much angular momentum as several smaller 
and more distant subhalos would, resulting in a decreased 
value of A s . 

The correlation between A m and f s is a significant re- 
sult. As A s decreases with increasing f s , the overall cor- 



relation between A and f s (see Figure 1240 must be en- 
tirely due to the increase in A m . We are therefore seeing 
evidence of the transfer of angular momentum from cap- 
tured subhalos, inducing a 'spin-up' of the mother (we 
re-emphasize that the mother is that largest structure 
in each halo, consisting of all the particles that do not 
belong to any substructures). There are two possible 
processes by which this can occur. In the first scenario, 
a recently absorbed dwarf halo is slowed by gravitational 
dynamical friction as it penetrates into the cluster. This 
results in a transfer of energy and angular momentum 
to the surrounding particles, whils t slow i ng the captured 
halo into a more stable orbit (see IZWl ll200l for a de- 
tailed study of this process). Hence, even mother halos 
that host a single massive subhalo near the halo cen- 
ter will have a greater spin due to the transfered an- 
gular momentum as t he su bhalo migrates towards the 
center. iPeirani et all {2004) and iVitvitska etafl l)2002fi 
directly measure this effect by comparing the radial an- 
gular momentum distri bution before and a fter a merger 
episode (see Figure 6 in IPeirani et alJ l)2004|) '). In general 
they find that the maximum specific angular momen- 
tum increases sharply in the afterm ath of th e merg er. 
This process is also seen indirectly bv lGardnerj (|2001^ in 
his simulations — halos that experienced a recent merger 
have a larger spin than those that have grown purely 
through gradual accretion. 

The second process by which substructure can transfer 
angular momentum to the host halo is through the tidal 
stripping of their outer regions as they orbit around the 
halo core. As Figure [23 shows, particles belonging to 
substructure have a greater angular momentum per unit 
mass than those that belong to the host halo. Therefore, 
as particles are stripped from substructure and become 
part of the mother halo, the mean angular momentum of 
the mother increases. 

6.4. Alignment of host and subhalo orbital angular 
momentum 
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In the preceding Section we demonstrated that clusters 
with a larger fraction of their mass in substructure have 
a greater spin. This is due to the higher average orbital 
angular momentum per unit mass of subhalos compared 
to the mother halo. However, for this to be true the 
majority of subhalos must have a prograde orbit with re- 
spect to the mother halo particles, i.e. the angular mo- 
mentum vector for the substructure must be well aligned 
with that of the host. In Figure TT7\ we investigate this 
by plotting the distribution of the relative orientations 
of the angular momentum vectors of substructure parti- 
cles (J s ), mother halo particles (3 m) and the cluster as 
a whole (3 to t — 3 s + Jm) for each halo in our sample. 
The middle panel (3 S ■ 3m) is of the most interest. It 
clearly demonstrates that there is a tendency for 3 S and 
3m to be closely aligned. Hence, it appears that sub- 
halos do tend to orbit in the same direction as the host 
halo particles. 

It has been demonstrate d both through obser- 
vations llPlionis fc Basil akps] | 2002|) and in sim- 
ulations (IColberg et all 120051 [Gill et alJ [2004b; 
iFaltenbacher et alJl2002j) that subhalos preferentially fall 
into cluste rs along t he filament ary structures surr ound- 
ing them. iTormenl l|1997l) and iGill et all (|2004bD both 
conclude that infall along filaments account for both the 
cluster shape and velocity structure. The correlation 
between the orientation of the angular momenta of 
mother halo and substructure particles is just a conse- 
quence of the same process responsible for the spinning 
up of a collapsing halo. The angular momentum of 
a halo is acquired through the influence of the tidal 
fields of neighboring structures acting on the infalling 
particles and subhal os during halo formation (see, for 
example Whitg [l984L and references therein), the influx 
of which is primarily from the high density regions that 
eventually form filamentary structures. Once captured, 
the outer region of a subhalo is tidally stripped, and the 
dense remnant spirals into and is eventually consumed 
by the halo core. Thus, the mother accumulates all the 
orbital angular momentum of the subhalo. Over time, 
as subhalos are captured and absorbed, the the spin 
of the mother halo gradually increases in this manner. 
Hence, the measured alignment between the angular 
momentum of the mother halo and subhalo particles 
is merely due to the temporal auto-correlation of the 
angular momenta of matter accumulated by the cluster. 

Additionally, a smaller halo on a prograde orbit around 
a cluster will gravitationally perturb nearby particles 
within the virial radius, causing a slight extension o f 
the halo shape towards it ijTremaine &: Weinberd ll984'). 
This in turn increases the attraction of the cluster on 
the halo, thus improving the likelihood that it will be 
captured. Halos that pass the cluster in retrograde or- 
bits will not be able to perturb the cluster to as great an 
effect. Hence, there is a slightly greater capture proba- 
bility for a halo that has a lower relative velocity to the 
local cluster particles, whether it arrives from a filament 
or not. 

6.5. Alignment of angular momentum and halo axes 

For a rotating body, it is expected that the direction 
of the angular momentum vector J will align with one 
of the principle axis vectors. The most stable configu- 
ration is when J is parallel to the minor axis c. Fig- 
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Fig. 27. — Distributions of the relative orientation of the angular 
momentum vectors for the entire cluster and the substructure (top 
panel), the substructure and the mother halo (middle panel) and 
the mother halo and entire cluster (bottom panel) . The middle plot 
demonstrates that there is reasonably good alignment (to within 
25°) of the angular momentum of the mother halo and subhalos. 



urc 1281 displays the distribution of the cosine of the an- 
gle 4> between the angular momentum vector and the 
major, intermediate and minor axes for all the halos in 
our sample. There is clearly a tendency for the angular 
momentum vector to be aligned with the minor axis of 
the halos. Correspondingly, the major axis is far more 
likely to be nearly perpendicular to J. This result was 
also fo und bvlBarnes fc Efsta thioul ill 9 871). I Warren et all 
lll992D.lBaTlin fc Steinmetzl (I2004D . lAlleood et alJ l)2005fl 
and lAvila-Reese et alJ l)2005|) .~ 

The intermediate axis alignment is interesting: it ap- 
pears to show no preferred alignment with J, having an 
almost uniform distribution over all cos0. We have in- 
vestigated whether this may be due to the difficulties in 
determining the minor axis when c/b is close to unity. 
However, we found this not to be the case — the axis 
ratios are widely distributed for these halos. Hence, a 
minority of halos are aligned so that their intermediate 
axes is parallel to the direction of the angular momentum 
vector — an unstable configuration for a solid rotating 
body. In Figure we break down the range of A into 
four equally sized sets. The purpose of this is to investi- 
gate whether the oddly aligned halos are those with the 
least (normalized) angular momentum. It is clear that 
as A decreases the distribution of orientations becomes 
more evenly spread. Therefore, the orientation of a halo 
is therefore dependant on the magnitude of its spin. The 
greater the angular momentum of a halo, the more likely 
it is to be found rotating around the minor principle axis, 
and the less likely around the major axis. Except for 
the most rapidly rotating halos, which maintain the ex- 
pected orientation, the results for the intermediate axis 
are inconclusive. The seemingly random distribution of 
alignments is maintained for slowly rotating halos. 

6.6. Correlation between A and triaxiality 

We now consider whether the IPeeblesI l|1969|) relation 
between the ratio of the major and minor principle axes 
and the cluster spin (Equation lll|) holds for the dark 
matter halos in our simulation. As we determined in 
Section 14.41 most of the halos tend to be prolate rather 
than oblate. Here we investigate whether there is an 
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cos(<|>) 

Fig. 28. — Distribution of the cosine of the angle <f> between 
the angular momentum vector and the major (solid), intermediate 
(dotted) and minor (dashed) axes. The plot is normalized so that 
a uniform distribution would have value of unity on the vertical 
axis (grayed horizontal line). 
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Fig. 29. — The distribution of the cosine of the angle between the 
angular momentum vector and the major (top panel), intermediate 
(middle), and minor (bottom) axis for halos grouped into ranges 
of total halo spin, A. The plot demonstrates an improvement in 
alignment as A increases. Rapidly spinning halos have a tendency 
for J to be more closely aligned with the minor axis. Each plot 
is normalized so that a uniform distribution would have value of 
unity on the vertical axis (grayed horizontal lines). 



increased departure from a spherical morphology (i.e. a 
reduction of the axis ratio c/a) as A increases. 

There are at least two effects that could influence this 
correlation. The first of these is the possibility of a sig- 
nificant misalignment between the minor principle axis 
and the angular momentum vector. Eauation ll ll assumes 
that the object is rotating around its minor principle 
axis. In those halos for which there is a significant mis- 
alignment between the minor axis and the angular mo- 
mentum vector, the effect of halo spin on shape will be 
unpredictable. However, as many of the halos in our sam- 
ple are prolate, with a mean c/b axis ratio near 0.9, we 
include only those halos for which the angle between the 
major axis and the angular momentum vector is greater 
than or equal to 70°. A second possible source of un- 
certainty will be the effects of substructure on the halo 
shape. Figure I23 showed the principle axis ratios using 
particles belonging only to the mother halo (i.e. having 



removed the subhalos) finding it to be more spherical 
than the entire halo. Equation II II applies to a cylmdri- 
cally symmetric and uniformly rotating polytrope, and 
therefore will certainly not hold when there is signifi- 
cant substructure in the cluster halo. Hence we discard 
all particles belonging to substructure before calculating 
the c/a axis ratio for each halo, and we use A m — the spin 
calculated using the mother halo particles only — rather 
than A. 

Even having accounted for these possible sources of 
discrepancy, we find that there is no correlation between 
the spin and the overall shape of the halos. It is there- 
fore clear that other processes are responsible for deter- 
mining halo morphology. The mean value we obtain 
for Ao is 0.033. Therefore, cluster mass halos cannot 
support themselves through rotation alone, it is instead 
through their velocity dispersion that they are able to 
do so. Hence, spin is largely irrelevant in determining 
shape even though we find a strong correlation between 
the angular momentum vector and the minor principle 
axis. We have also shown that the angular momentum 
of dark matter halos is increased as subhalos are accu- 
mulated, which also influences halo shape. Furthermore, 
subhalos tend to be found in prograde orbits, partly be- 
cause of a preferential direction of infall from filaments 
surround the halo. It appears that the orientation of dark 
matter clusters can be explained by the specific merging 
directions of smaller halos that they capture, during their 
formative years. 

7. DISCUSSION AND CONCLUSION 

In this study, we set out to investigate the overall phys- 
ical characteristics of a sample of over 2000 massive dark 
matter cluster halos identified in a 320 h- Mpc, 1024 3 
particle simulation . Using the algorithm described in 
iWeller et all l)2005|) . all particles belonging to substruc- 
ture in each cluster were identified, down to a minimum 
subhalo mass of 30 particles, or ~ 10 11 /i~ 1 Af Q . We then 
applied two further criteria to our sample. First, to en- 
sure that all the halos in our dataset were well resolved — 
thus ensuring that no substructure has been washed out 
by the so-called 'over-merging' problem — we set a min- 
imum cluster halo mass of 10,000 particles, leaving us 
with an overall sample of 2159 halos. Second, we identi- 
fied and tagged those halos not yet in a state of dynamical 
equilibrium. The aim of the second criterion is to enable 
separate analyses of those halos that are virialized and 
those that are still at an earlier stage of their formation. 

In order to implement this second criterion, we defined 
a measure of virialization, (3 (Equation El, based on the 
virial theorem and including a necessary term to com- 
pensate for the surface pressure at the virial radius of 
each halo. Using this technique, we demonstrated quan- 
titavely that higher mass, substructure-rich halos are less 
virialized. We removed from our sample any halo with 
< —0.2 (i.e. those still collapsing) - which make up 
roughly 3.4% of the entire sample - tagging them as 'un- 
virialized'. In total, our sample contained 2085 virialized 
halos, each containing between 10,000 and 500,000 parti- 
cles. Such a large number of well resolved halos allowed 
us to determine their physical properties in a statistically 
significant way. 

We then constructed a catalogue of their physical prop- 
erties, including concentration, morphology, angular mo- 
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mentum, maximum circular velocity, and the fraction 
of mass in substructures. The main questions motivat- 
ing our analysis were: how are the halo properties dis- 
tributed, both in our main halo sample and in those that 
we have tagged as 'unvirialized' ? Do halo properties de- 
pend on mass? And how are they influenced by substruc- 
ture? These questions can only be answered if one anal- 
yses a sample of considerable size with no prior selection 
of specific halos to resimiulate from a lower resolution 
simulation. The distributions obtained in this work are 
drawn from a large enough sample (^2000) to address 
these questions in a statistical meaningful way. 

The halos in our virialized sample contained a median 
value of 5.6% of their mass in substructure, with a tail to 
large substructure fractions. On average, they tended to 
have a prolate morphology and a mean minor to major 
principle axis ratio of c/a = 0.707 ± 0.095. The unviri- 
alized halos were more substructure-rich, with a median 
value of 19.6% of their mass in substructure. They also 
tended to be less spherical, have a higher spin and a 
lower concentration than their dynamically relaxed coun- 
terparts. 

We have clearly demonstrated that the physical prop- 
erties of halos are dependant on their mass and therefore 
their age. The median fraction of mass contained in sub- 
structure increases with halo mass as f s cx (M v ; r /-M*) m , 
where m = 0.44 ± 0.06. Th is supports the prediction 
made bv lZentner et al.l l)2005f) that as higher mass halos 
have formed more recently than lower mass halos, they 
should contain more substructure as they have had less 
time to disrupt their internal subhalo population. One 
direct consequence of this is that higher mass halos ap- 
pear more prolate: we find c/a cx (M vir /M*) s , where 
5 = —0.049 ± 0.007. As expected, halo concentration de- 
creases as mass increases, with c(M) oc (M/Af») Q with a 
= -0.12 ± 0.03 , agree ing well with the s l opes m easured 
bv lDolag et all IpTOlJ) and iBullock et~aD l|2001bD . 

Halos with greater amounts of substructure tend also 
to have a higher spin; the angular momentum per unit 
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mass of subhalos is far greater than that of the mother 
halo. We also find evidence of the transfer of angular 
momentum to the mother halo (which consists of all the 
particles in the halo that do not belong to substructures) 
from captured subhalos, via dynamical friction and the 
tidal stripping of their outer regions. Furthermore, the 
orbital angular momentum of subhalos is typically well 
aligned with that of their host. We propose that this is 
mainly due to the process by which a halo accumulates 
angular momentum from captured subhalos, which, due 
to large-scale tidal fields, have a preferential direction of 
infall from the surr ounding filaments (see, for example 
IColberg et alJ 12005(1 . The angular momentum of con- 
sumed subhalos is added to that of the mother halo, 
so that over time it will begin to rotate in the same 
direction to which the subhalos are 'injected', explain- 
ing the observed alignment. Additionally, we find that 
the orientation of the angular momentum vector of halo 
with respect its morphology is dependant its magnitidue. 
Rapidly spinning halos tend to align with the minor prin- 
ciple axis; halos with very low values of A are have a 
higher probabiliy of being aligned with the intermedi- 
ate or major axes. However, angular momentum is not 
responsible for determining halo shape. Indeed, we con- 
cluded that the correlation between angular momentum 
and the minor principle axis is due to a common influence 
on both properties — the preferential merging direction 
of subhalos along halo filaments. 
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